clear
load('IOdata.mat');
n=size(X',1);m=size(X,1);s=size(Y,1);
epsilon=10^-10;
f=[zeros(1,n) -epsilon*ones(1,m+s) 1];
A=zeros(1,n+m+s+1);b=0;
LB=zeros(n+m+s+1,1);UB=[];
LB(n+m+s+1)=-inf;
for i=1:n
   Aeq=[X eye(m) zeros(m,s) -X(:,i)
        Y zeros(s,m) -eye(s) zeros(s,1)];
    beq=[zeros(m,1)
        Y(:,i)];
    w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);  
end
w;
lambda=w(1:n,:);
s_minus=w(n+1:n+m,:);
s_plus=w(n+m+1:n+m+s,:);
theta=w(n+m+s+1,:);


